Freezing of an unconventional two-dimensional plasma 
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We study an unconventional two-dimensional, two-component classical plasma on a sphere, with 
emphasis on detecting signatures of melting transitions. This system is relevant to Ising-type quan- 
tum Hall states, and is unconventional in the sense that it features particles interacting via two dif- 
ferent two-dimensional Coulomb interactions. One species of particles in the plasma carries charge 
of both types (Qi,Q2), while the other species carries only charge of the second type (0, — Qi)- 
We find signatures of a freezing transition at Q\ ~ 140. This means that the species with charge 
of both types will form a Wigner crystal, whereas the species with charge of the second type also 
shows signatures of being a Wigner crystal, due to the attractive inter-component interaction of the 
second type. Moreover, there is also a Berezinskii-Kosterlitz-Thouless phase transition at Q\ ~ 4, 
at which the two species of particles bind to form molecules that are neutral with respect to the 
second Coulomb interaction. These two transitions appear to be independent of each other, giving a 
rectangular phase diagram. As a special case, Q2 — describes the (conventional) two-dimensional 
one-component plasma. Our study is consistent with previous studies of this plasma, and sheds new 
light on the freezing transition of this system. 

PACS numbers: 73.43. Cd, 74.20.De, 74.25.Uv 



I. INTRODUCTION 

Multi-component quantum condensates with novel 
types of inter-component interactions are of consider- 
able interest in contemporary physics. For example, 
they are relevant to widely disparate systems, including 
low-dimensional spin-1/2 quantum antifcrromagnets 1 5 , 
Bose- Einstein condensates^"^, multi-component /multi- 
band superconductors^^] and non-Abelian quantum 
Hall states and topological superconductors^. These 
systems have the remarkable property of possessing a 
mapping to a classical multi-component plasma system 
with highly unusual intra- and inter-component inter- 
actions. The statistical properties of these unconven- 
tional plasmas - especially their phase diagrams - have 
important ramifications for the physics of their corre- 
sponding fractional quantum Hall systems^. The sta- 
tistical physics of such systems has only recently begun 
to be explored. In a previous paper, we investigated the 
metal-insulator transition in a particular version of such 
a plasmE^. In this paper, we will extend these investi- 
gations to a study of the freezing of such a plasma from 
a liquid to a Wigner crystal. 

The canonical partition function of the unconven- 
tional two-component plasma that we investigate is given 

bySsnu 

Z = j (ll d ^) (n d a «.) «" V . « 



where the potential energy 

N N 

V= -Ql ln|w a -w b | + Ql ln|zi-w„| 

a<b— 1 a, z—l 

N 

- (Ql + Ql) E ln l z * + ^-bg ( 2 ) 
i<j=l 

describes two species (components) of particles inter- 
acting via two different types of two-dimensional (2D) 
Coulomb interactions, which are logarithmic. Here, the 
Zj are coordinate vectors for the N particles of component 
z, which carry charge Q\ of the first interaction (type 1) 
and charge Q2 of the second interaction (type 2). The w a 
are coordinate vectors for the iV particles of component 
w, which carry no charge of type 1 and charge —Q2 of 
type 2. The term V z ,bg describes the interaction of the z- 
particles with a uniform density neutralizing background 
charge. Note that the form given in Eq. implies that 
the temperature T = 1. 

This plasma is related^ to inner products of quantum- 
mechanical trial wave functions of Ising-type quan- 
tum Hall s tates, such as the Moore-Read PfaffiarfSl, 
anti-Pfaffia n 1 15 * 16 !, and Bonderson-Slingerland hierarchy 
states^. For the charge values relevant to these states, 
the plasma was shown to be in its metallic liquid phase^, 
which allows for the calculation of the braiding statistics 
of quasiparticle excitations of these states 1 -^, confirming 
their conjectured non-Abelian statistics. This plasma 
is also related to rotating two-component Bose-Einstein 
condensates (BECs) in two dimensions^. 
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While Ref. 13 focused on the cases Q\ — 0, 2, which are 
particularly relevant for Ising-type quantum Hall states, 
here we will investigate the plasma for large values of Q ± . 
In the limit in which Q2 = 0, the w-particles do not in- 
teract, and the plasma thus reduces to the standard 2D 
one-component Coulomb plasma (OCP). It is generally 
believed that, at high values of Q\, the OCP will be in 
a 2D solid state in which the charges form a triangu- 
lar lattice with quasi-long-range translational and long- 
range orientational order^, as found in the simulations 
in Refs. IT5H2"2l However, some studies have claimed that 
there is no low-temperature (high-Qi) crystalline state 
in the OCP, due to the proliferation of screened disclina- 
tionpHH. 

If we assume that the generally-held view is correct 
(and we present evidence supporting this view), so that 
there is a low-temperature crystalline state, then the 
melting of this crystal can occur according to either of 
two possible scenarios. One possibility is the Kosterlitz- 
Thouless-Halperin-Nelson- Young theory (KTHNYjPHS], 
according to which dislocation pairs unbind via a vector- 
defects version of the BKT transition at a Berezinskii- 
Kosterlitz-Thouless (BKT) like transition. The system 
then enters a hexatic liquid phase in which there is no 
translational order, but there is quasi-long-range hexatic 
order. Then, there is a second BKT transition at which 
disclination pairs unbind, hexatic order is lost, and the 
system enters an isotropic liquid phase. The other possi- 
bility is a direct first-order melting transition at a lower 
temperature than the KTHNY-theory predicts^. There 
have been considerable efforts to investigate 2D melting, 
both experimentally and by numerical simulations. Some 
studies have found KTHNY transitions while others have 
found a weakly first-order melting transitiorP^SI. It ap- 
pears that the nature of 2D melting depends on details of 
the interatomic potential. In the case of logarithmic in- 
teractions, most numerical simulations find a first-order 
transitiorpSHUl 

Before proceeding to a description of our simulations, 
we mention that, in principle, there is one other possibil- 
ity: a Lifshitz transition from the liquid to a striped or 
"microemulsion" phase and then later to a Wigner crys- 
tal, as discussed by Kivelson and SpivalPO. Such a sce- 
nario must be considered when there is a linear coupling 
between the order parameter and the uniform density 
(i.e., between the order parameter at wavevector q and 
the density at wavevector — q) or, equivalently, when the 
first derivative of the energy with respect to the density 
is discontinuous at the transition. However, in our case, 
the order parameter is the density at non-zero wavevec- 
tor, so no such linear coupling can occur. Furthermore, 
the order parameter vanishes on both sides of the tran- 
sition since the crystalline phase is only quasi-long-range 
ordered, so there would be no discontinuity even if there 
were a linear coupling. However, even in systems to 
which the Kivelson- SpivalP^ argument applies, there are 
two possible scenarios, similar to the ones that we con- 
sider: a direct first-order phase transition (which is per- 



mitted for the case of logarithmic interactions) and a con- 
tinuous transition via one or more intermediate phases. 



II. MODEL AND SIMULATION 

The system described in Eqs. ([!]) and ^ is studied by 
means of large-scale Monte Carlo simulations on a sphere 
of radius R. In this geometry, the distance between two 
points ri and r 2 is taken to be the chord length 



| ri -r 2 | = V2#(l-fi-f 2 )*, 



(3) 



and the term V z ,bg is simply a uniform constant that 
can be disregarded. Hence, the model in Eq. ([2|) may be 
written in the form (up to constant terms ^flll 2 ^ 40 ! 41 ! 
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- (Qi + Ql) E ln d 

i<j=l 



N 



E Mi 
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Zj • Zj) 



W a • Wj,) 

(4) 



Here, w a , Zj are the positions of the particles on the sur- 
face of the unit sphere. Details of the derivation, as well 
as on the technicalities of the Monte Carlo simulations, 
are presented in Ref. 1131 Moreover, to improve sampling 
at high value s of Q±, we used the parallel tempering al- 
gorithrrp214£y w j iere the set of couplings was found by 
measuring first-passage-times as described in Ref. l44l 

In addition to the logarithmic interactions, we regu- 
larize the attractive interactions by adding a short-range 
hard-core repulsion such that particles are not permitted 
to be closer than the particle diameter d. Hence, there 
is a nonzero dimensionless density 77 = 2Ns/A where 
s = ircP/A and A is the area of the system. 



III. RESULTS FOR THE ONE-COMPONENT 
PLASMA 

First, we consider the case in which Q2 = 0. This is 
motivated by the fact that previous studies of the OCP 
on the surface of a sphere are not consistent. In Ref.[20l a 
freezing transition at Q\ ~ 140 was found by comparing 
the free energy of the solid and liquid state. However, 
in Ref. [UJ the absence of a finite-temperature crystalline 
state was claimed and numerical evidence supporting this 
was provided, essentially by showing that the correlation 
length for crystalline order was non-divergent, £ cx Q\ 
for all Qi 1 > 0. We will return to this below. 

The structure function is given by 



1 

N 
= U 



dr dr'e iq -( r - r ')(n(r)n(r')) 



N ^ 



q-(ri-r 3 -)\ 



(5) 
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Figure 1: (Color online) Correlation length £ as a function of 
Qi for 8 different system sizes in the range 100 < N < 4800. 
The inset is a plot of S(q) for the specific case when N = 3200 
and Qi = y/UO. 

Here, (n(r)n(r')) is the density-density function, where 
n(r) — J2j 5{r — Tj) is the local density, and (...) de- 
notes a statistical average. Furthermore, {r^} denotes 
the positions of the point particles in the problem and 
q is the Fourier space vector. In this work, we measure 
the azimuthal ave rage of 5 (q) modified for a spherical 
geometry, given bySSESHSl 

S(q) = 1 + 2vmi? 2 / dd [g{R9) - 1] sin 9J (qR6), (6) 
Jo 

where n is the number density, R is the radius of the 
sphere on which the particles live, g(R9) is the pair dis- 
tribution function with 9 as the chord angle, Jq(x) is a 
zeroth order Bessel function, and q is the magnitude of 

q- 

The inset of Fig. [I] shows a plot of S(q). We assume 
that the correlation length £ is inversely proportional to 
the width of the first peak in S(q), and may thus be de- 
termined by a Lorentzian fit. The procedure is identical 
to that used in Ref. [351 and the result is given in Fig. [T] 
For small values of Q\ our results are similar to Fig. 2 
in Ref. 25, However, when Qi ps 12, a value that cor- 
responds well with the critical coupling of the freezing 
transition, we find a kink developing with increasing N, 
that clearly violates £ cx Q±. This kink is not seen in 
Fig. 2 of Ref. [25j However, we note that the markers 
of that figure exhibits large scattering. Moreover, the 
authors did not consider larger values of Q±. 

A hallmark of a 2D solid is that translational correla- 
tions have a power-law decay, (e lG ( r-r )) ~ |r — r'|~' ,G , 
where G is the reciprocal lattice vector, r, r' arc lattice 
points in the 2 D solid and t]g is a temperature dependent 
exponent! 29 * 46 ^. Consequently, the first-order Bragg peak 
in 5(q) will, as seen by Eq. (|5]), scale as 5(G) ~ L 2 ~ ,,G , 
where L cx A 1 / 2 is the spatial linear extent of the sys- 
tem. Now, by integrating over the Bragg peak of a 2D 



solicP^, the finite-size scaling of the azimuthally averaged 
first peak in S(q), is given by 

5(G) ~ L 1 -^ ~ at(i-xg)/2. ( 7 ) 

Fig.[2]shows the results for the maximum value of the first 
peak in S(q) for a wide range of system sizes and for dif- 
ferent values of Q\. As for £, we find that the peak value 
also exhibits a kink at Q\ ps 140 that should be associ- 
ated with an abrupt change in the translational correla- 
tions in the plasma. Indeed, when studying the finite-size 
behavior more closely in the lower panel of Fig. [2j the re- 
sults show that when Q\ < 130, 5(G) ~ const, when 
TV increases. This is the behavior expected in the liquid 
phase, with exponentially decaying translational correla- 
tions where 5(G) ~ £ 2 . However, when Q\ > 150, the 
results clearly show that there is a positive slope that 
develops with increasing A, thus confirming the finite- 
size behavior of the 2D solid given in Eq. 0. When 
Q\ = 140, it is difficult to determine whether the system 
is in the solid phase or not, suggesting that Q\ = 140 is 
close to the melting point of the OCP. Note that in Fig. [2] 
the height of the first-order peak in 5(g), 5(G) ps 5 when 
Q\ = 140. This is consistent with the 2D freezing cri- 
terion for a crystal with long-range int eract ions (charac- 
terized by a divergent bulk modulu a 18 ^ 20 * 48 *. 

A key prediction of the KTHNY theory is that rj^ < 
1/3 in the solid phase, where the limiting value of 1/3 is 
reached at the critical point of melt ing from a triangular 
lattice to the hexatic phas#^2129| j± s a resu ^ j m this 
scenario, 5(G) grows more rapidly with N than A 1 / 3 for 
all Q\ greater than the critical value; 5(G) grows as N 1 ^ 3 
at the transition point; and 5(G) saturates in the liquid 
phase. Meanwhile, if the transition were first-order, the 
limiting value of ?7g would be smaller than 1/3, so that 
5(G) would grow more rapidly than N 1 / 3 at the tran- 
sition point, i.e., the slowest possible growth of 5(G) in 
the crystalline phase would be faster than N 1 ^ 3 . Conse- 
quently, we expect the slope of In 5(G) vs. In AT to be 
steeper than 1/3 for all Q\ in the crystalline phase or, 
by the results above, for all Q\ > 140. By determining 
the slope of In 5(G) vs. In A at Q\ ps 140, we could then 
determine if the transition is of KTHNY type or is first- 
order. However, as may be seen in Fig. [2] the slopes 
of In 5(G) vs. In A in the putative crystalline phase 
arc not steeper than 1/3 in our simulations. However, 
the slopes steepen with increasing A, possibly converg- 
ing towards the expected behavior in the thermodynamic 
limit. Therefore, we are unable to determine which type 
of transition occurs, nor whether melting proceeds via an 
intermediate hexatic phase (see below). 

It is worth emphasizing that for Q\ < 140, 5(G) ap- 
pears to saturate to a finite value, as expected in a liquid, 
while, for Q\ > 140, 5(G) does not appear to saturate, as 
expected in a crystal (although, as noted above, it does 
not grow as rapidly as expected). Therefore, the lower 
panel of Fig.[2]is also qualitatively consistent with a crys- 
talline phase of the 2D OCP, which melts at Q\ ps 140. 
Taken together with Fig. [T] and the upper panel of Fig. [2] 
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Figure 2: (Color online) Finite-size behavior of the peak value 
of the structure factor, S(G), as a function of coupling Ql and 
size N. The upper panel shows S(G) as a function of Ql for 
eight different sizes in the range 100 < N < 4800. The lower 
panel is a log-log plot of S(G) as a function of N for seven 
fixed values of Ql. The dashed line is a reference line that 
yields the expected finite-size behavior at the melting point 
according to KTHNY theory. Lines are guide to the eyes. 



orientational susceptibility 



X6 =n 




Here, we have 



*6 = 



1 



N 

_ V e i6( ^ 



Furthermore, the quantity 



1 '»« 



(8) 



(9) 



(10) 



may be thought of as the local bond-orientational order 
parameter of particle i, where all bond angles 4>i a are 
measured with respect to the closest nearest neighbor in 
the tangential plane of particle i. The sum is over all rii 
nearest neighbors as determined by Voronoi construction. 
In Eq. |9]), the angle is the bond angle of the bond to 
the closest nearest neighbor of particle i, measured with 
respect to the i,j chord in the tangential plane of parti- 
cle i. Thus, the chord line combining the two particles 
for every term in the sum of Eq. |9]) , serves as a line of 
reference for bond-orientational ordeiP^^l 

In Fig.[3j the results for X6 & re given. When the system 
size is large (N > 800), a peak is found. However, the 
value of Qi at which the peak occurs, appears to converge 
towards Ql « 140. Thus, with the resolution available, 
we cannot confirm the existence of a hexatic phase since 
the position of the peak in X6 does not appear to converge 
to a coupling significantly different from Q\ « 140. Our 
findings are consistent with earlier works, which also have 

found no traces of an intermediate hexatic phase in these 
systemsSsHll 



IV. RESULTS FOR AN UNCONVENTIONAL 
TWO-COMPONENT PLASMA 



this provides clear evidence for the existence of a low- 
temperature crystalline state of the OCP on a sphere, in 
agreement with previous studiesSsHUl. Our results con- 
tradict the claims made in Ref . [5S] for the non-existence 
of a crystalline phase. The present work considers larger 
values of Q\ and provides considerably better statistics 
than Ref. 012 

The transition from the hexatic to the isotropic liq- 
uid phase is governed by fluctuations in the orientational 
order parameter. These can be computed by the bond- 
orientational susceptibility which exhibits a peak at this 
transition. When this peak is obtained at a different 
point than the onset of quasi-long-range translational or- 
der, it indicates the existence of an intermediate hex- 
atic phas e in b etween the crystalline and isotropic liq- 
uid phase^ES! Qn a sphere, we compute the bond- 



We now turn our attention to the full model in Eq. ([2| , 
i.e., when both Q\ and Q2 are nonzero. In particular, we 
consider how the translational ordering of both w and 
z particles is affected as we increase the coupling con- 
stant in the second interaction-channel, Qi. As for the 
OCP, we study translational correlations by measuring 
the structure factors S w (q), S z (q) defined by Eq. (|6| with 
S(q), g(R0) — !• S w / Z (q), g w / z {R0). In addition, we also 
measure the inverse dielectric constant for charges with 
interaction of type 2, given by 



= 1 - 



irQ 2 2 R 2 
A 




(11) 



This quantity measures the screening properties for 
charges interacting with Q2, and it signals a charge- 
unbinding transition involving z- and w-particles^. 
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Figure 3: (Color online) Bond-orientational susceptibility \6 
as a function of Qi for 7 different system sizes in the range 
100 < N < 3200. The position of the peak appears to con- 



verge to Qi w 140. 



In Fig. |4j results are given for the height of the first- 
order peak in the structure factor for component z and 
w, for the case when Q\ = \. Apart from the fact that 
the height of the peak in the structure factor is much 
larger for the z particles than the w-particles, the size- 
and Qf-dependence of the peaks are qualitatively very 
similar for the two components. In particular, they both 
exhibit a kink at Q\ s» 140, which should be associated 
with melting of a 2D solid, similar to the OCP case in 
the upper panel of Fig. [2] Specifically, when we extract 
the finite-size behavior in the log-log plots in Fig. [5j we 
find that both components exhibit S(G) ~ const., con- 
sistent with being in the liquid phase, when Q\ < 130. 
When Qi > 150, the results clearly show that there is a 
power-law dependence on N, consistent with the finite- 
size behavior of a 2D solid. These results are consistent 
with the phase diagram in Fig. [6j 

The inverse dielectric constant e^ 1 is measured to be 
zero to the left of the red line in Fig. [6} Thus, the w and z 
particles are in a metallic state regardless of the change 
in the structural properties when Q\ m 140. This is 
most salient with respect to the second type of Coulomb 
interaction (which has an effective strength that is de- 
termined by e^)- I n the liquid phase, it is clear that w 
and z particles are in a metallic state. In the crystalline 
phase, there are interstitials and vacancies in the crystal, 
so that a finite fraction of w and z particles should be 
considered as unbound particles that are able to screen 
test particles interacting with charges of type 2, thereby 
leading to e^ 1 =0. At larger values of there is 
a transition at which w and z particles are bound into 
molecules. For Q\ above this transition point, which is at 
Q\ c ~ 4, €22 has a non-zero value, as found in Ref. [131 
Although the w particles are able to screen the type 2 
interaction when z particles form a Wigner crystal, their 
translational correlations exhibit signatures of a 2D solid 
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Figure 4: (Color online) Results from the Monte Carlo sim- 



ulations when Q2 = 1 and 77 



10- 



Panel (a) shows the 



height of the first-order peak in the structure factor for the 
z particles, S Z (G), as a function of the coupling Q\ for seven 
different sizes in the range 100 < N < 3200. Panel (b) shows 
the height of the first-order peak in the structure factor for 
the w particles, S W (G), as a function of the coupling Q\ for 
seven different sizes in the range 100 < N < 3200. In order 
to give an impression on how a typical structure factor looks 
like, the insets of panel (a) and (b) show plots of S z (q) and 
Sw{q) for the specific case when N = 800, Q\ = 180 and 
Ql = 1. 



(as seen in Fig [4]), attributed to a higher probability of 
the w particles to be co-centered with z particles due 
to the attractive inter-component interactions of type 2. 
On average, a finite fraction of the w particles should 
be considered as bound to the z particles, thus adapt- 
ing to the 2D crystalline structure that is created by the 
strong repulsive interactions among the z particles, when 
Qi > 140. The signatures of freezing of the w-particles 
is thus an effect which is induced by the freezing of the 
z-particles. 

Hence, we can summarize the situation as follows, as 
depicted in Fig. [6j In Phase I, the w and z particles are 
unbound and are separately in a liquid state. In Phase II, 
the w and z particles are bound into molecules that are 
neutral with respect to the second type of Coulomb inter- 
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Figure 5: (Color online) Log-log plots of the results in panel 
(a) and (b) in Fig. [4] Both panels show the height of the 
first-order peak in the structure factor as a function of size 
N. In the lower panel, S' W (G) is the height of the peak of 
S m (G) when we have subtracted the regular part in order to 
properly extract the singular finite-size behavior of S w (q) in a 
log-log plot. The solid lines are reference lines that yield the 
expected finite-size behavior at the melting point according 
to KTHNY theory. Lines are guide to the eyes. 



action, and these molecules form a liquid. In Phase III, 
the w and z particles are unbound; the z particles form a 
Wigner crystal while the w particles form a liquid, albeit 
one with modulated density due to its interaction with 
the Wigner crystal. In Phase IV, the w and z particles 
are bound into molecules forming a Wigner crystal. In 
Appendix[5] we explain the details of how the transition 
lines were obtained. 

We consider phase III to be a 2D counterpart of the 
situation that was reported for a three-dimensional sys- 
tem in Ref. [8] This work considered a two-component 
rotating BEC with a negative dissipationless Andrccv- 
Bashkin drag^. It was found that in this mixture, a 
situation may arise where the component with the small- 
est stiffness will be a modulated vortex liquid. That is, the 
soft component breaks translational symmetry while ex- 
hibiting an unbroken symmetry in order parameter space. 
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Figure 6: (Color online) The phase diagram as a function 
of Ql and Q%. The dashed red line is the charge unbinding 
transition, at which z and w particles become bound together 
for at Q\ above this line. The dotted blue line is the melting 
line of the Wigner crystal. In Phases I and II, the z particles 
are in a liquid state; in Phases III and IV, the z particles 
form a Wigner crystal. In Phases I and III, the w and z 
particles are unbound; in Phases II and IV, they are bound 
into molecules comprised of one z and one w particle. See the 
text for details. 



The vortices of the soft component are likely to be co- 
centered with the vortices of the stiffest component, and 
will thus adapt to the spatial structure of the latter. As 
shown in Ref. [T31 the MR plasma corresponds to a 2D 
two-component rotating BEC with negative drag, where 
the z component is stiffer than the w component when 
Ql > 0. 



V. SUMMARY AND CONCLUSIONS 

In summary, we have considered the melting of an 
unconventional 2D two-component plasma on a sphere 
with particles interacting in two different channels, which 
may be viewed as an analogous plasma describing a non- 
Abelian Ising-type quantum Hall state or a realization 
of a two-component two-dimensional Bose-Einstein con- 
densate with inter-component non-dissipative drag. In 
the limiting case where there are no interactions of type 
2 (Q2 — 0), the system is a standard 2D one- component 
plasma. Both for the one-component plasma and the 
unconventional two-component plasma, we find that the 
system freezes on a sphere for large enough inter-particle 
interactions. For the two-component plasma, the w- 
particles do not have strong intra-component interac- 
tions, but still show signatures of forming a 2D solid. 
This is attributed to the attractive inter-component in- 
teractions with the z-particles that leads to a higher 
probability of the w-particles to be co-centered with z- 
particles. The w-particles nonetheless form a metallic 
state. We have also examined the possible existence of an 



7 



intermediate hexatic phase in the one-component plasma. 
Our results show that the value of Q\ where diclinations 
unbind and orientational order is lost, cannot be distin- 
guished from the value of Qi where dislocations unbind 
and quasi-long-range translational order is lost. 
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Appendix A: Determination of the transition lines 

We now discuss the determination of the phase transi- 
tion lines in Fig. [6] in more detail. 

First, consider the line (red in Fig. [6]) at which the z 
and w particles unbind. In this work, we find the critical 
point of the BKT transition by curve-fitting the inverse 
dielectric constant to a logarithmic finite-size scaling rela- 
tion with one free parameter (see Appendix C in Ref . IT3")) . 
This means that we assume that the transition is a BKT 
transition as we use the BKT value of the universal jump 
in the finite-size scaling relation. Thus, it is a slightly less 
self-consistent approach than what was used in Ref. [T3J 
but still, one can regard this as a verification of the BKT 
nature, as one should not expect a good fit to the scaling 
relation if the transition is of a different nature.^ 

The curve-fitting was performed according to the de- 
scription in Appendix C in Ref. [TBI for sizes N = 70, 100, 
150, 200, 300, and 500, for Q\ = 20, 100, and 160 and 
for densities 77 = 0.001, 0.0004, and 0.0001. In Fig. [7) 
the results for the transition point Q\ c as a function of r\ 
and Qi are given. We have also included the results for 
Q\ = and 2 from Ref. 13 In order to obtain a crude 
estimate of the transition temperature in the low density 
limit, we extrapolate to 77 = by fitting the results for 



finite i] to a power law Q\ (rf) 



ari b where Q\ , 



4.4 



c? 4 - 2 



Qi = -+ 

QI = 2 — 
Q? = 20 

Qj = 100 * 

Q\ = 160 h~° 



0.0005 
'/ 



0.001 



Figure 7: (Color online) The transition point is determined 
for rj — > by extrapolating the Q\ c values obtained for non- 
zero 77, assuming a power-law dependence. Lines are guide to 
the eyes. 



a, and b are free parameters. The estimates we find are: 



QL 


= 4.016 ±0.002 


for 


Qi 


= 


01,0 


= 4.015 ± 0.004 


for 


Qi 


= 2 


QL 


= 4.013 ±0.101 


for 


Qi 


= 20 


QL 


= 4.012 ±0.060 


for 


Qi 


= 100 


Qlc 


= 3.963 ±0.070 


for 


Qi 


= 160 



(Al) 



These values are plotted in the phase diagram in Fig. |6j 
and we take the phase boundary to be the best fit straight 
line running through them. 

We now consider the Wigner crystal melting transi- 
tion, depicted by the blue line in Fig. [6] This transition 
is found by measuring the value of Qi a t which S z (G) 
attains its maximum second-derivative. (See Fig. [2] and 
Fig. [4] in the paper for example.) In Fig. [8j we show 
the estimates of the transition point Q\ function 
of inverse system size TV -1 for Q\ = 0, 1, 3, and 5. The 
transition points are estimated by averaging the results 
for N > 800, with errors determined by a bootstrap anal- 
ysis. The estimates we find are: 



Q 2 lc = 140.6 ± 1.5 
Qf' c = 140.3 ±2.0 
Qf' c = 144.4 ±4.9 
Ql' c = 142.6 ±8.2 



for Ql = 

for Ql = 1 

for Ql = 3 

for Ql = 5 



(A2) 



The phase diagram in Fig. [6] is obtained by using these 
values, and we take the phase boundary to be the best 
fit straight line running through them. 
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Figure 8: (Color online) The transition point Q\ c as a func- 
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